Equations (3) to (5) describe the INS error dynamics using a set of nonhomogeneous differential equations with time-varying coefficients. The full determination of the INS error propagation is too complicated when taking into account the real trajectories and maneuvers. In this case, we generally use a simulation approach to aid the analysis (Titterton and Weston , 2004; Groves , 2008). The situation becomes much more complex for an aided INS because the integrated navigation solution is affected not only by INS error sources but also by the external measurement noises. Fortunately, in the railway track surveying application, the nominal trajectory and motion of the host track trolley are simple and almost deterministic, which reduces the complexity and makes the analytic assessment possible.
Figure 2 depicts a flowchart of the semi-analytic assessment of the relative measurement accuracy with the aided INS in measuring railway track irregularity. The procedure consists of an analytic assessment phase (the upper part of the figure) and a Monte Carlo simulation phase (the lower part). First, the system model and measurement model are simplified according to the several assumptions that hold in the specific railway surveying scenario. Then the coupling between the channels of the INS error model vanishes, and the error state dynamics is reduced to a linear time-invariant stochastic system driven only by white Gaussian noise, as enclosed by the dashed lines in the figure. The simplified INS error differential equations are solved in the Laplace domain, where \({\mathbf {H}}_\text {INS}(s)\) represents the transfer function. The Laplace transformation of the steady-state KF estimator is performed to yield the optimal estimate, where \({\mathbf {H}}_\text {KF}(s)\) characterizes the transfer function. An analytical solution of the navigation errors in the Laplace domain, i.e., \(\delta \hat{\varvec{r}}(s)\), is obtained as the difference between \(\delta \varvec{r}(s)\) and \(\hat{\delta \varvec{r}}(s)\).
The Monte Carlo simulation part, enclosed in the gray box, is performed as follows: create a transfer function model based on the derived \(\delta \hat{\varvec{r}}(s)\) model; then, use the simulated white noise series as the input to generate the navigation error samples, i.e., the system responses; and finally, evaluate the relative measurement accuracy based on the output navigation error samples. More details on the simulation part are given in subsequent sections.
Fig. 2Flowchart of the semi-analytical analysis of the relative measurement accuracy of the GNSS/INS integrated system
Full size imageSimplification for the railway measurement case Assumption 1The track trolley is assumed to move uniformly in a straight line in the south-north direction.
This assumption is based on the following facts: High-speed railway track is usually designed with a very large radius of curvature in both the horizontal and vertical profiles and with a slope gradient smaller than 1.5%. The longest chord in railway track irregularity measurements does not exceed 300 m; thus, it is reasonable to simplify such a short curve section to a straight line and to assume that the host vehicle moves uniformly in south-north direction and remains level, in which case the attitude matrix \({\mathbf {C}}_b^n = {\mathbf {I}}\). As defined in appendix, the horizontal railway track deformation refers to the deviation of rails from its nominal position in the lateral and vertical directions. When the trolley moves in the south-north direction, we can evaluate the horizontal and vertical irregularity measurement accuracy by directly analyzing the east and vertical positioning errors, respectively. Thus, this assumption would simplify the following analysis. In addition, the trolley moves at a low speed, and its attitude changes slowly, which makes the dynamics-induced errors (such as the effects of the scale factor and cross-coupling) negligible.
Assumption 2The IMU is mounted with its sensitive axes perfectly aligned with the host vehicle axes, i.e., the b-frame is coincident with the v-frame. In this case \({\mathbf {C}}_b^v\) becomes the identity matrix.
This assumption is reasonable because the IMU mounting angles can be estimated and compensated for with sufficient accuracy, as discussed in (Chen et al. , 2020)
Assumption 3Lever arms of the GNSS antenna and odometer are all zeros, and GNSS positions with centimeter accuracy, obtained by the postprocessed kinematic, are available all the time, as the lever arms can be measured with sufficient accuracy and the related effect can be corrected in practice.
Assumption 4The local gravity uncertainty is assumed to be much smaller than the accelerometer measurement error and is therefore negligible.
System model simplificationThe performance of the aided INS is also influenced by the real trajectories and maneuvers. We list the key information on the trajectory, inertial sensor errors, and absolute navigation error as follows:
mean position of the trajectory: latitude = 30\(^\circ\), longitude = 114\(^\circ\), h = 20 m.
mean absolute positioning error: \(\delta r_N\) = 0.01 m, \(\delta r_E\) = 0.01 m, \(\delta r_D\) = 0.02 m.
velocity: \(v_N\) = 1 m/s, \(v_E\) = 0 m/s, \(v_D\) = 0 m/s.
mean velocity error: \(\delta v_N\) = 0.002 m/s, \(\delta v_E\) = 0.002 m/s, \(\delta v_D\) = 0.002 m/s.
mean attitude error: \(\phi _N\) = 0.003\(^\circ\), \(\phi _E\) = 0.003\(^\circ\), \(\phi _D\) = 0.005\(^\circ\).
gyro bias: \(\delta \omega _{ib}^b\) = 0.01 \(^\circ /h\).
accelerometer bias: \(\delta f^b\) = 10 mGal.
local gravity value: 9.78 m/s\(^2\).
We can calculate the magnitude of each term in the INS error differential equations by substituting the parameters listed above. The error terms with a magnitude smaller than 10% of the predominant term are considered insignificant terms and thus are neglected in the following analysis. The inertial sensor error terms \(\delta \varvec{\omega }_{ib}^b\) and \(\delta \varvec{f}^b\) are always kept. As a result, the INS error differential equations from (3) to (5) can be simplified as
$$\begin{aligned} {\left\{ \begin{array}{ll} \delta {\dot{r}}_N = \delta v_N \\ \delta {\dot{r}}_E = \delta v_E \\ \delta {\dot{r}}_D = \delta v_D \\ \delta {\dot{v}}_N = -f_D \phi _E + \delta f_N \\ \delta {\dot{v}}_E = f_D \phi _N + \delta f_E \\ \delta {\dot{v}}_D = \delta f_D \\ {\dot{\phi }}_N = -\delta \omega _{ib,N}^n \\ {\dot{\phi }}_E = -\delta \omega _{ib,E}^n \\ {\dot{\phi }}_D = -\delta \omega _{ib,D}^n \end{array}\right. } \end{aligned}$$(23)where \(\delta f_N\), \(\delta f_E\), \(\delta f_D\) are the accelerometer measurement errors in the north, east, and vertical directions, respectively. \(\delta \omega _{ib,N}^n\), \(\delta \omega _{ib,E}^n\), \(\delta \omega _{ib,D}^n\) are the gyroscopic measurement errors in the north, east, and vertical directions, respectively. After simplification, as shown in (23), the coupling between channels vanishes, and each channel can be analyzed separately.
Measurement model simplificationA similar simplification can be carried out for measurement Eqs. (10) and (14) by taking into account the assumptions and trajectory information listed above, yielding
$$\begin{aligned} \varvec{z}_r = \delta \varvec{r}^n + \varvec{n}_{r} \end{aligned}$$(24)$$\begin{aligned} \varvec{z}_\text {v} = \delta \varvec{v}^n - \varvec{v}^n \times \varvec{\phi } + \varvec{n}_{v} \end{aligned}$$(25)Measurement error propagation modelThe navigation error propagation modeling of the GNSS/INS is analyzed in the vertical and horizontal channels. In the following, the details on the vertical channel analysis are presented, while the analysis of the horizontal channel is similar and given in the appendix.
For the vertical channel analysis, in addition to the height and vertical velocity errors, the easting gyro bias \(\delta b_{gE}\) and vertical accelerometer bias \(\delta b_{aD}\) should be augmented into the error state vector. Additionally, since the NHC is induced as a navigation aid, the attitude error term \(\phi _E\) should be augmented. Thus, the error state vector and the related state dynamics model can be written as
$$\begin{aligned} \varvec{x}_D= & {} \left[ \begin{matrix} \delta h&\delta v_D&\phi _E&\delta b_{gE}&\delta b_{aD} \end{matrix}\right] ^\text {T} \end{aligned}$$(26)$$\begin{aligned} \dot{\varvec{x}}_D(t)= & {} {\mathbf {F}}_D(t)\varvec{x}_D(t) + {\mathbf {G}}_D(t)\varvec{w}_D(t) \end{aligned}$$(27)with
$$\begin{aligned} \begin{matrix} {\mathbf {F}}_D=\left[ \begin{matrix} 0&{} -1&{} 0&{} 0&{} 0\\ 0&{} 0&{} 0&{} 0&{} 1\\ 0&{} 0&{} 0&{} -1&{} 0\\ 0&{} 0&{} 0&{} -1/T_{gb}&{} 0\\ 0&{} 0&{} 0&{} -1/T_{ab}&{} 0\\ \end{matrix} \right] \text {, }&{} {\mathbf {G}}_D=\left[ \begin{matrix} 0&{} 0\\ 0&{} {\mathbf {I}}_4\\ \end{matrix} \right] \\ \end{matrix} \end{aligned}$$(28)where the subscript ‘D’ denotes the down direction and \(\varvec{w}_D(t)\) represents the driving white Gaussian noise of strength \({\mathbf {Q}}_D\):
$$\begin{aligned} \varvec{w}_D= & {} \left[ \begin{matrix} 0&w_{aD}&w_{gE}&w_{gbE}&w_{abD} \end{matrix}\right] ^\text {T} \end{aligned}$$(29)$$\begin{aligned} E\big \{\varvec{w}_D(t) \varvec{w}^\text {T}_D(t+\tau ) \big \}= & {} {\mathbf {Q}}_D \delta (\tau ) \end{aligned}$$(30)where \(\delta b_{gE}\) is modeled as a first-order Gauss–Markov process with correlation time \(T_{gb}\) and driven by white Gaussian noise process \(w_{gbE}\). The residual bias in the vertical accelerometer \(\delta b_{aD}\) is modeled as a first-order Gauss–Markov process with correlation time \(T_{ab}\) and driven by white Gaussian noise process \(w_{abD}\). \(w_{aD}\) is the white Gaussian noise that corrupted the vertical accelerometer measurement. \(w_{gE}\) is the white Gaussian noise that corrupted the gyro measurement in the east direction.
The measurements in the KF for the aided INS in the vertical channel includes the height difference between the INS and GNSS and the vertical velocity difference between the INS and the NHC. According to (25), the velocity measurement equation can be written as
$$\begin{aligned} z_{vD} = \delta v_D + v_E \phi _N- v_N \phi _E + n_{vD} \end{aligned}$$(31)where \(n_{vD}\) is the vertical component of \(\varvec{n}_{v}\) in (12) and represents the measurement uncertainty of the across-track zero velocity. According to assumption 1, the host vehicle is assumed to move in the south-north direction, in which case the east velocity shall be zero. Thus, (31) can be simplified as
$$\begin{aligned} z_{v,D} = \delta v_D - v_N \phi _E + n_{vD} \end{aligned}$$(32)Therefore, the measurement equation for the aided INS in the vertical channel can be written as
$$\begin{aligned} \varvec{z}_D = {\mathbf {H}}_D \varvec{x}_D + \varvec{n}_D \end{aligned}$$(33)where \({\mathbf {z}}_D = \left[ \begin{matrix} z_h&z_{vD}\end{matrix}\right] ^\text {T}\), \({\mathbf {n}}_D = \left[ \begin{matrix} n_{rD}&n_{vD}\end{matrix}\right] ^\text {T}\)
$$\begin{aligned} {\mathbf {H}}_D= & {} \left[ \begin{matrix} 1 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} -v_N &{} 0 &{} 0 \end{matrix}\right] \end{aligned}$$(34)$$\begin{aligned} E\big \{\varvec{n}_D(t) \varvec{n}^\text {T}_D(t+\tau )\big \}= & {} {\mathbf {R}}_D \delta (\tau ) \end{aligned}$$(35)However, if the vehicle moves with a nonzero constant heading angle rather in the south-north direction, then \(v_E \ne 0\). In this case the simplification from (31) to (32) is not valid and the measurement Eq. (31) should be used.
The state dynamics equations of the continuous KF are given by
$$\begin{aligned} \dot{\hat{\varvec{x}}}_D(t) = {\mathbf {F}}_D(t) \hat{\varvec{x}}_D(t) + {\mathbf {K}}_D(t) \left[ \varvec{z}_D(t) - {\mathbf {H}}_D(t) {\hat{x}}_D (t)\right] \end{aligned}$$(36)It is obvious from Eqs. (28) and (34) that \({\mathbf {F}}_D\), \({\mathbf {G}}_D\), \({\mathbf {H}}_D\) are all constant matrices; \({\mathbf {R}}_D\) and \({\mathbf {Q}}_D\) are constant and diagonal matrices. The filter gains converge to steady-state values in a short time, resulting in constant matrices \({\mathbf {P}}\) and \({\mathbf {K}}\). The matrix \({\mathbf {P}}\) can be calculated by solving the continuous matrix Riccati differential equation \(\dot{{\mathbf {P}}}(t) = {\mathbf {0}}\). The gain matrix \({\mathbf {K}}_D\) in this case is a 5-by-2 matrix.
Taking the Laplace transform of Eq. (27) yields the related solution in the Laplace domain:
$$\begin{aligned} \varvec{x}_D(s) = \left[ \left( s{\mathbf {I}} - {\mathbf {F}}_D\right) ^{-1} {\mathbf {G}}_D \varvec{w}_D(s) \right] \end{aligned}$$(37)The expansion of the above formula can be written as
$$\begin{aligned} \delta h(s)= & {} -\frac{1}{s^2} w_{aD}(s) - \frac{1}{s^2(s+\gamma _{ab})} w_{abD}(s) \end{aligned}$$(38a)$$\begin{aligned} \delta v_D(s)= & {} \frac{1}{s} w_{aD}(s) + \frac{1}{s (s+\gamma _{ab})}w_{abD}(s) \end{aligned}$$(38b)$$\begin{aligned} \phi _E(s)= & {} \frac{1}{s} w_{gE}(s) - \frac{1}{s (s+\gamma _{gb})} w_{gbE}(s) \end{aligned}$$(38c)$$\begin{aligned} \delta b_{aD}(s)= & {} \frac{1}{s+\gamma _{ab}} w_{abD}(s) \end{aligned}$$(38d)$$\begin{aligned} \delta b_{gE}(s)= & {} \frac{1}{s+\gamma _{gb}} w_{gbE}(s) \end{aligned}$$(38e)where \(\gamma _{gb} = \frac{1}{T_{gb}}\), \(\gamma _{ab} = \frac{1}{T_{ab}}\). By taking the Laplace transform of Eq. (33) and substituting (37) into the transformed equation we obtain \(\varvec{z}_D(s)\), whose expansion can be written as
$$\begin{aligned} z_h(s)= & {} -\frac{1}{s^2} w_{aD}(s) - \frac{1}{s^2(s+\gamma _{ab})} w_{abD}(s) \nonumber \\&+ n_{rD}(s) \end{aligned}$$(39a)$$\begin{aligned} z_{vD}(s)= & {} \frac{1}{s} w_{aD}(s) - \frac{v_N}{s} w_{gE}(s) + \frac{1}{s(s+\gamma _{ab})} w_{abD}(s) \nonumber \\&+ \frac{v_N}{s(s+\gamma _{gb})} w_{gbE}(s) + n_{vD}(s) \end{aligned}$$(39b)Solving Eq. (36) in the Laplace domain yields
$$\begin{aligned} \hat{\varvec{x}}_D(s) = \left[ \left( s{\mathbf {I}} - {\mathbf {F}}_D + {\mathbf {K}}_D {\mathbf {H}}_D\right) ^{-1} {\mathbf {K}}_D \right] \varvec{z}_D(s) \end{aligned}$$(40)Subtracting the height component \(\hat{\delta h}(s)\) in vector \(\hat{\varvec{x}}_D(s)\) from \(\delta h(s)\) in (37) yields the height estimate error of the aided INS, \(\delta {\hat{h}}(s)\), which is a linear function of \(\varvec{w}_D(s)\), written as
$$\begin{aligned} \begin{aligned} \delta {\hat{h}}(s) = ~&H_{aD}~w_{aD}(s) + H_{gE}~ w_{gE}(s) \\&+ H_{gbE} ~w_{gbE}(s) + H_{abD}~ w_{abD}(s) \\&+ H_{nrD}~ n_{rD}(s) + H_{nvD} ~n_{vD}(s) \end{aligned} \end{aligned}$$(41)where \(H_{aD}\), \(H_{gE}\), \(H_{gbE}\), \(H_{abD}\), \(H_{nrD}\), and \(H_{nvD}\) are the related coefficients of the noise terms. The above equation is the error propagation model of the aided INS in the vertical channel in the Laplace domain. It tells that the height error of the GNSS/INS integration for measuring the railway track irregularity can be regarded as a time-invariant stochastic system with the white noise as the system input. The detailed expression of \(\delta {\hat{h}}(s)\) is a bit complicated, and the related derivation can be got by a symbolic operation with the MATLAB software.
Relative measurement accuracy analysisThe position error solution \(\delta {\hat{r}}_\text {est}(s)\) in the Laplace domain has been obtained. Theoretically, the corresponding analytic solutions in the time domain, denoted by \(\delta {\hat{r}}_\text {est}(t)\), can be obtained by an inverse Laplace transform. Then, we can completely evaluate the relative measurement accuracy of the aided INS according to (A.5) and (A.7). Unfortunately, obtaining this analytic expression \(\delta {\hat{r}}_\text {est}(t)\) is impossible because the inverse Laplace transformation of white noise in Eqs. (41) and (B.12) does not lead to an explicit analytic expression. Therefore, the remainder analysis should be performed with a Monte Carlo simulation to generate the samples of \(\delta {\hat{r}}_\text {est}(t)\), as depicted in Fig. 2.
The Monte Carlo simulation relies on repeated random sampling to obtain numerical results, which is helpful in understanding the behaviors of a stochastic system that are not amenable to analysis by the usual direct mathematical methods (Brown and Hwang , 2012). Since white Gaussian noise is a stationary and Gaussian-distributed random process with a constant spectral density function over all frequencies, if it is put into a linear time-invariant system, the system output will also be stationary and Gaussian distributed. It is evident from (41), (B.12) and (20) that \(\triangledown \hat{\varvec{r}}_\text {est}\) will also be a stationary, Gaussian-distributed, and ergodic random process. In that case, we would be able to conduct the performance analysis of the relative measurement accuracy with sufficiently long samples \(\delta \hat{\varvec{r}}_\text {est}(t)\). The analysis procedure is as follows:
(1)Use the MATLAB function tf to create a transfer function model based on the error propagation model (41) and (B.12) in the Laplace domain.
(2)Simulate the input white noise samples, including the system noise, measurement noise, and driving noise of the first-order Gauss–Markov process; the corresponding parameters needed are listed subsequently.
(3)Use the MATLAB function lsim to generate the time response of the dynamic system (generated in step 1) to the input stimuli (the output in step 2), obtaining navigation error samples of the aided INS in the time domain.
(4)Calculate the railway track irregularity measurement error \(\triangledown \hat{\varvec{r}}_\text {est}\) in (20) and its corresponding covariance matrix \({\mathbf {P}}_\text {ir}\) in (21).
The information needed for generating the input white noise are listed as follows:
Angular Random Walk (ARW) = 0.002 \(^\circ /\sqrt{h}\), Velocity Random Walk (VRW) = 0.001 \(m/s/\sqrt{h}\).
\(T_{gb}\) = 1000 s, \(\sigma _{gb}\) = 0.005 \(^\circ /h\), \(T_{gb}\) = 1000 s, \(\sigma _{gb}\) = 25 mGal
\(n_{rE}\) = 1 cm/\(\sqrt{Hz}\), \(n_{rD}\) = 2 cm/\(\sqrt{Hz}\), \(n_{v}\) = 0.1 mm/s/\(\sqrt{Hz}\)
Figure 3 is a sample of the aided INS positioning error process in the lateral and vertical directions. In this figure, the navigation errors are plotted versus the travel distance, which is obtained from the calibrated odometer sensor. Since the positioning error in the along-track direction has little influence on the track irregularity measurement accuracy, it is not included in this figure. The first observation is that the error sequence varies within 1 cm and shows an explicit spatial correlation.
Based on the generated navigation error samples, we can calculate the railway track irregularity measurement error \(\triangledown \hat{\varvec{r}}_\text {est}\) as introduced in the above step 4. Figure 4 is a histogram of \(\triangledown \hat{\varvec{r}}_\text {est}\) in the east and vertical directions, which are related to the short-wave alignment and longitudinal direction, respectively, illustrating the distribution of \(\triangledown \hat{\varvec{r}}_\text {est}\). They approximately follow the Gaussian distribution, as expected. The results show that 1 mm accuracy is achievable in vertical and horizontal track irregularity measurements with a navigation-grade IMU aided by a carrier-phase differential GNSS and the NHC. This result answers the first question raised in the introduction section.
Fig. 3A sample of the aided INS positioning error process in the lateral and vertical directions
Full size imageFig. 4Histogram of the short-wave track irregularity measurement error incurred by the aided INS
Full size image